## This code include scripts for Appendix Figures
# Figure A3-a and A3-b
# Figure A4-a and A4-b
# Figure A5-a, A5-c, and A5-b
# Figure A6-a, A6-c, and A6-b

rm(list = ls())

library(readr)
library(ggplot2)
library(stargazer)
library(lme4)
library(dplyr)
library(haven)
library(texreg)
library(tidyr)

load("Final_JPR/data/paperdata.RData")
data <- paperdata %>%
       filter(tek_count>0)

source("Final_JPR/Code/coefplot.R")
source("Final_JPR/Code/marginal_effects.R")
source("Final_JPR/Code/setup.R")




############################################ Figure A3a-A3b
## drop population less than 500k

data_500km <- data %>% dplyr::filter(pop_total_linear >500000)
dv <- 'onset_do_flag' 
## controls
ctrs <- c('lineq_total', 'status_excl','lsize','family_warhistdummy','family_downgraded2', 
          'ctygdppc_log', 'ctypop_log','peaceyears','peaceyears_sq','peaceyears_cub')

ivs_1 <- c('best_tlineq_total', 'tekbest_status_excl')
ivs_2 <- c('worst_tlineq_total', 'tekworst_status_excl')
ivs_3 <- c("median_tlineq_total",'tekmedian_status_excl')
ivs_4 <- c('nearst_tlineq_total','teknearst_status_excl')

f_do_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))

#
ivs <- c('best_tlineq_total','worst_tlineq_total','median_tlineq_total','nearst_tlineq_total')


## using log of transnational inequality
model_onset_500k <- list()
for (i in 1:length(ivs)){
       model_onset_500k[[i]] <- glm(eval(parse(text = paste0('f_do_', i))), 
                                    data = data_500km, family = binomial(link = "logit"))
}

p1 <- marginaleffect_logit(ModelResults = model_onset_500k, n.sim = 1000, 
                           varname = c("best_tlineq_total", "worst_tlineq_total",
                                       "median_tlineq_total","nearst_tlineq_total"),
                           data = data_500km,
                           val1 = 0, val2 =1, clusterid = "gwgroupid") 

p1 +  scale_x_continuous(breaks = 1:4,
                         labels = parse(text =rev(c( 
                                "Nearest~TEK",
                                "Median~TEK",
                                "Worst~TEK",
                                "Best~TEK")))) + ggtitle("Conflict Onset")

ggsave("Final_JPR/Figures/appendix_figures/fig_A3-a.pdf", width = 6.5, height = 4.2, units='in', dpi=600)

############################################ Figure A3b
## dependent variable
dv <- 'incidence_flag' 
## controls
ctrs <- c('lineq_total', 'status_excl','lsize','family_warhistdummy','family_downgraded2', 
          'ctygdppc_log', 'ctypop_log','peaceyears','peaceyears_sq','peaceyears_cub')

ivs_1 <- c('best_tlineq_total', 'tekbest_status_excl')
ivs_2 <- c('worst_tlineq_total', 'tekworst_status_excl')
ivs_3 <- c("median_tlineq_total",'tekmedian_status_excl')
ivs_4 <- c('nearst_tlineq_total','teknearst_status_excl')

f_do_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))

#
ivs <- c('best_tlineq_total','worst_tlineq_total','median_tlineq_total','nearst_tlineq_total')


## using log of transnational inequality
model_incidence_500k <- list()
for (i in 1:length(ivs)){
       model_incidence_500k[[i]] <- glm(eval(parse(text = paste0('f_do_', i))), 
                                        data = data_500km, family = binomial(link = "logit"))
}



p2 <- marginaleffect_logit(ModelResults = model_incidence_500k, n.sim = 1000, 
                           varname = c("best_tlineq_total", "worst_tlineq_total",
                                       "median_tlineq_total","nearst_tlineq_total"),
                           data = data_500km,
                           val1 = 0, val2 =1, clusterid = "gwgroupid") 

p2 + scale_x_continuous(breaks = 1:4,
                        labels = parse(text =rev(c( 
                               "Nearest~TEK",
                               "Median~TEK",
                               "Worst~TEK",
                               "Best~TEK"))))+ ggtitle("Conflict Incidence")

ggsave("Final_JPR/Figures/appendix_figures/fig_A3-b.pdf", width = 6.5, height = 4.2, units='in', dpi=600)



############################################################ Figure A4a-A4b

### relative terms
dv <- 'onset_do_flag' 
## controls
ctrs <- c('lineq_total', 'status_excl','lsize','family_warhistdummy','family_downgraded2', 
          'ctygdppc_log', 'ctypop_log','peaceyears','peaceyears_sq','peaceyears_cub')

ivs_1 <- c('best_low_tlineq', 'best_high_tlineq', 'tekbest_status_excl')
ivs_2 <- c('worst_low_tlineq', 'worst_high_tlineq', 'tekworst_status_excl')
ivs_3 <- c('median_low_tlineq', 'median_high_tlineq','tekmedian_status_excl')
ivs_4 <- c('nearst_low_tlineq', 'nearst_high_tlineq','teknearst_status_excl')

f_do_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))

#
ivs <- c('best','worst','median','nearst')

## using log of transnational inequality
Relmodel_onset_500k <- list()
for (i in 1:length(ivs)){
       Relmodel_onset_500k[[i]] <- glm(eval(parse(text = paste0('f_do_', i))), 
                                       data = data_500km, family = binomial(link = "logit"))
}


p2_1 <- marginaleffect_logit2(ModelResults =  Relmodel_onset_500k, n.sim = 1000, 
                              varname = c("best_low_tlineq", "worst_low_tlineq",
                                          "median_low_tlineq","nearst_low_tlineq"),
                              data = data_500km,
                              dv = "relative_poverty",
                              val1 = 0, val2 =1, clusterid = "gwgroupid") 


p2_2 <- marginaleffect_logit2(ModelResults =   Relmodel_onset_500k, n.sim = 1000, 
                              varname = c("best_high_tlineq", "worst_high_tlineq",
                                          "median_high_tlineq","nearst_high_tlineq"),
                              data = data_500km,
                              dv = "relative_wealth",
                              val1 = 0, val2 =1, clusterid = "gwgroupid") 


df_relative_onset <- bind_rows(p2_1, p2_2)



p = ggplot(df_relative_onset , aes(colour = dv)) + 
       geom_hline(yintercept = 0, lty = 2) +
       geom_linerange(aes(x = id, ymin = low,
                          ymax = high),
                      lwd = 1.5, position = position_dodge(width = 1/2)) +
       geom_pointrange(aes(x = id, y =median, ymin =low,
                           ymax =high, shape = dv),
                       lwd = 1.5, position = position_dodge(width = 1/2),
                       fill = "WHITE")+
       geom_text(aes(x = id +0.1 ,
                     y = median,
                     label = format(round(median, 2))),position = position_dodge(width = 1/2)) +
       theme_bw() + xlab("") + ylab("First Difference in Predicted Prob.") + 
       theme(legend.position="bottom",
             legend.title=element_blank(),
             legend.text = element_text(colour="blue", size=12),
             axis.text= element_text(size=12),
             text = element_text(size=12),
             plot.title = element_text(hjust = .5, size = 12))+
       ggtitle("Conflict Onset")+
       scale_colour_discrete(labels = rev(c("Relative Wealth","Relative poverty")))+
       scale_shape_discrete(labels = rev(c("Relative Wealth","Relative poverty")) )+
       guides(color = guide_legend(nrow = 1))+
       scale_x_continuous(breaks = 1:4,labels =parse(text =rev(c( 
              "Nearest~TEK",
              "Median~TEK",
              "Worst~TEK",
              "Best~TEK")))) 

ggsave("Final_JPR/Figures/appendix_figures/fig_A4-a.pdf", width = 6.5, height = 4.2, units='in', dpi=600)


## dependent variable
dv <- 'incidence_flag' 
## controls
ctrs <- c('lineq_total', 'status_excl','lsize','family_warhistdummy','family_downgraded2', 
          'ctygdppc_log', 'ctypop_log','peaceyears','peaceyears_sq','peaceyears_cub')

ivs_1 <- c('best_low_tlineq', 'best_high_tlineq', 'tekbest_status_excl')
ivs_2 <- c('worst_low_tlineq', 'worst_high_tlineq', 'tekworst_status_excl')
ivs_3 <- c('median_low_tlineq', 'median_high_tlineq','tekmedian_status_excl')
ivs_4 <- c('nearst_low_tlineq', 'nearst_high_tlineq','teknearst_status_excl')

f_do_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))

#
ivs <- c('best','worst','median','nearst')

## using log of transnational inequality
Relmodel_incidence_500k <- list()
for (i in 1:length(ivs)){
       Relmodel_incidence_500k[[i]] <- glm(eval(parse(text = paste0('f_do_', i))), 
                                           data = data_500km, family = binomial(link = "logit"))
}


p1_1 <- marginaleffect_logit2(ModelResults = Relmodel_incidence_500k, n.sim = 1000, 
                              varname = c("best_low_tlineq", "worst_low_tlineq",
                                          "median_low_tlineq","nearst_low_tlineq"),
                              data =data_500km,
                              dv = "relative_poverty",
                              val1 = 0, val2 =1, clusterid = "gwgroupid") 


p1_2 <- marginaleffect_logit2(ModelResults =  Relmodel_incidence_500k, n.sim = 1000, 
                              varname = c("best_high_tlineq", "worst_high_tlineq",
                                          "median_high_tlineq","nearst_high_tlineq"),
                              data =data_500km,
                              dv = "relative_wealth",
                              val1 = 0, val2 =1, clusterid = "gwgroupid") 


df_relative_incidence <- bind_rows(p1_1, p1_2)



p = ggplot(df_relative_incidence , aes(colour = dv)) + 
       geom_hline(yintercept = 0, lty = 2) +
       geom_linerange(aes(x = id, ymin = low,
                          ymax = high),
                      lwd = 1.5, position = position_dodge(width = 1/2)) +
       geom_pointrange(aes(x = id, y =median, ymin =low,
                           ymax =high, shape = dv),
                       lwd = 1.5, position = position_dodge(width = 1/2),
                       fill = "WHITE")+
       geom_text(aes(x = id +0.1 ,
                     y = median,
                     label = format(round(median, 2))),position = position_dodge(width = 1/2)) +
       theme_bw() + xlab("") + ylab("First Difference in Predicted Prob.") + 
       theme(legend.position="bottom",
             legend.title=element_blank(),
             legend.text = element_text(colour="blue", size=12),
             axis.text= element_text(size=12),
             text = element_text(size=12),
             plot.title = element_text(hjust = .5, size = 12))+
       ggtitle("Conflict Incidence")+
       scale_colour_discrete(labels = rev(c("Relative Wealth","Relative poverty")))+
       scale_shape_discrete(labels = rev(c("Relative Wealth","Relative poverty")) )+
       guides(color = guide_legend(nrow = 1))+
       scale_x_continuous(breaks = 1:4,labels =parse(text =rev(c( 
              "Nearest~TEK",
              "Median~TEK",
              "Worst~TEK",
              "Best~TEK")))) 

ggsave("Final_JPR/Figures/appendix_figures/fig_A4-b.pdf", width = 6.5, height = 4.2, units='in', dpi=600)



############################################################ Figure A5a

### drop TEK with one
############################################ 
data_tek <- data %>% dplyr::filter(tek_count >1)
dv <- 'onset_do_flag' 
## controls
ctrs <- c('lineq_total', 'status_excl','lsize','family_warhistdummy','family_downgraded2', 
          'ctygdppc_log', 'ctypop_log', 'peaceyears','peaceyears_sq','peaceyears_cub')

ivs_1 <- c('best_tlineq_total', 'tekbest_status_excl')
ivs_2 <- c('worst_tlineq_total', 'tekworst_status_excl')
ivs_3 <- c("median_tlineq_total",'tekmedian_status_excl')
ivs_4 <- c('nearst_tlineq_total','teknearst_status_excl')

f_do_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))

#
ivs <- c('best_tlineq_total','worst_tlineq_total','median_tlineq_total','nearst_tlineq_total')


## using log of transnational inequality
model_onset_tek <- list()
for (i in 1:length(ivs)){
       model_onset_tek [[i]] <-glm(eval(parse(text = paste0('f_do_', i))), 
                                   data = data_tek, family = binomial(link = "logit"))
}




## dependent variable
dv <- 'incidence_flag' 
## controls
ctrs <- c('lineq_total', 'status_excl','lsize','family_warhistdummy','family_downgraded2', 
          'ctygdppc_log', 'ctypop_log', 'peaceyears','peaceyears_sq','peaceyears_cub')

ivs_1 <- c('best_tlineq_total', 'tekbest_status_excl')
ivs_2 <- c('worst_tlineq_total', 'tekworst_status_excl')
ivs_3 <- c("median_tlineq_total",'tekmedian_status_excl')
ivs_4 <- c('nearst_tlineq_total','teknearst_status_excl')

f_do_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))

#
ivs <- c('best_tlineq_total','worst_tlineq_total','median_tlineq_total','nearst_tlineq_total')


## using log of transnational inequality
model_incidence_tek <- list()
for (i in 1:length(ivs)){
       model_incidence_tek[[i]] <-glm(eval(parse(text = paste0('f_do_', i))), 
                                      data = data_tek, family = binomial(link = "logit"))
}


p1_1 <- marginaleffect_logit2(ModelResults =  model_onset_tek, n.sim = 1000, 
                              varname = c("best_tlineq_total", "worst_tlineq_total",
                                          "median_tlineq_total","nearst_tlineq_total"),
                              data = data_tek,
                              dv = "onset",
                              val1 = 0, val2 =1, clusterid = "gwgroupid") 


p1_2 <- marginaleffect_logit2(ModelResults =  model_incidence_tek, n.sim = 1000, 
                              varname = c("best_tlineq_total", "worst_tlineq_total",
                                          "median_tlineq_total","nearst_tlineq_total"),
                              data = data_tek,
                              dv = "incidence",
                              val1 = 0, val2 =1, clusterid = "gwgroupid") 


df_tek <- bind_rows(p1_1, p1_2)



p = ggplot(df_tek , aes(colour = dv)) + 
       geom_hline(yintercept = 0, lty = 2) +
       geom_linerange(aes(x = id, ymin = low,
                          ymax = high),
                      lwd = 1.5, position = position_dodge(width = 1/2)) +
       geom_pointrange(aes(x = id, y =median, ymin =low,
                           ymax =high, shape = dv),
                       lwd = 1.5, position = position_dodge(width = 1/2),
                       fill = "WHITE")+
       geom_text(aes(x = id +0.1 ,
                     y = median,
                     label = format(round(median, 2))),position = position_dodge(width = 1/2)) +
       theme_bw() + xlab("") + ylab("First Difference in Predicted Prob.") + 
       theme(legend.position="bottom",
             legend.title=element_blank(),
             legend.text = element_text(colour="blue", size=12),
             axis.text= element_text(size=12),
             text = element_text(size=12),
             plot.title = element_text(hjust = .5, size = 12))+
       ggtitle("Logistic Models (More than One TEK)")+
       scale_colour_discrete(labels = rev(c("Conflict Onset","Conflict Incidence")))+
       scale_shape_discrete(labels = rev(c("Conflict Onset","Conflict Incidence")) )+
       guides(color = guide_legend(nrow = 1))+
       scale_x_continuous(breaks = 1:4,labels =parse(text =rev(c( 
              "Nearest~TEK",
              "Median~TEK",
              "Worst~TEK",
              "Best~TEK")))) 

ggsave("Final_JPR/Figures/appendix_figures/fig_A5-a.pdf", width = 7.5, height = 4.2, units='in', dpi=600)


############################################################ Figure A5c
## dependent variable
dv <- 'onset_do_flag' 
## controls
ctrs <- c('lineq_total', 'status_excl','lsize','family_warhistdummy','family_downgraded2', 
          'ctygdppc_log', 'ctypop_log','peaceyears','peaceyears_sq','peaceyears_cub')

ivs_1 <- c('best_low_tlineq', 'best_high_tlineq', 'tekbest_status_excl')
ivs_2 <- c('worst_low_tlineq', 'worst_high_tlineq', 'tekworst_status_excl')
ivs_3 <- c('median_low_tlineq', 'median_high_tlineq','tekmedian_status_excl')
ivs_4 <- c('nearst_low_tlineq', 'nearst_high_tlineq','teknearst_status_excl')


f_do_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))

#
ivs <- c('best','worst','median','nearst')


## using log of transnational inequality
model_onset_relative_TEK <- list()
for (i in 1:length(ivs)){
       model_onset_relative_TEK[[i]] <- glm(eval(parse(text = paste0('f_do_', i))), 
                                            data = data_tek, family = binomial(link = "logit"))
}


## dependent variable
dv <- 'incidence_flag' 
## controls
ctrs <- c('lineq_total', 'status_excl','lsize','family_warhistdummy','family_downgraded2', 
          'ctygdppc_log', 'ctypop_log','peaceyears','peaceyears_sq','peaceyears_cub')

ivs_1 <- c('best_low_tlineq', 'best_high_tlineq', 'tekbest_status_excl')
ivs_2 <- c('worst_low_tlineq', 'worst_high_tlineq', 'tekworst_status_excl')
ivs_3 <- c('median_low_tlineq', 'median_high_tlineq','tekmedian_status_excl')
ivs_4 <- c('nearst_low_tlineq', 'nearst_high_tlineq','teknearst_status_excl')


f_do_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))

#
ivs <- c('best','worst','median','nearst')


## using log of transnational inequality
model_incidence_relative_TEK <- list()
for (i in 1:length(ivs)){
       model_incidence_relative_TEK[[i]] <- glm(eval(parse(text = paste0('f_do_', i))), 
                                                data = data_tek, family = binomial(link = "logit"))
}


p1_1 <- marginaleffect_logit2(ModelResults =  model_incidence_relative_TEK, n.sim = 1000, 
                              varname = c("best_low_tlineq", "worst_low_tlineq",
                                          "median_low_tlineq","nearst_low_tlineq"),
                              data = data_tek,
                              dv = "relative_poverty",
                              val1 = 0, val2 =1, clusterid = "gwgroupid") 


p1_2 <- marginaleffect_logit2(ModelResults =  model_incidence_relative_TEK, n.sim = 1000, 
                              varname = c("best_high_tlineq", "worst_high_tlineq",
                                          "median_high_tlineq","nearst_high_tlineq"),
                              data = data_tek,
                              dv = "relative_wealth",
                              val1 = 0, val2 =1, clusterid = "gwgroupid") 


df_relative_incidence <- bind_rows(p1_1, p1_2)



p = ggplot(df_relative_incidence , aes(colour = dv)) + 
       geom_hline(yintercept = 0, lty = 2) +
       geom_linerange(aes(x = id, ymin = low,
                          ymax = high),
                      lwd = 1.5, position = position_dodge(width = 1/2)) +
       geom_pointrange(aes(x = id, y =median, ymin =low,
                           ymax =high, shape = dv),
                       lwd = 1.5, position = position_dodge(width = 1/2),
                       fill = "WHITE")+
       geom_text(aes(x = id +0.1 ,
                     y = median,
                     label = format(round(median, 2))),position = position_dodge(width = 1/2)) +
       theme_bw() + xlab("") + ylab("First Difference in Predicted Prob.") + 
       theme(legend.position="bottom",
             legend.title=element_blank(),
             legend.text = element_text(colour="blue", size=12),
             axis.text= element_text(size=12),
             text = element_text(size=12),
             plot.title = element_text(hjust = .5, size = 12))+
       ggtitle("Conflict Incidence(More than One TEK)")+
       scale_colour_discrete(labels = rev(c("Relative Wealth","Relative poverty")))+
       scale_shape_discrete(labels = rev(c("Relative Wealth","Relative poverty")) )+
       guides(color = guide_legend(nrow = 1))+
       scale_x_continuous(breaks = 1:4,labels =parse(text =rev(c( 
              "Nearest~TEK",
              "Median~TEK",
              "Worst~TEK",
              "Best~TEK")))) 

ggsave("Final_JPR/Figures/appendix_figures/fig_A5-c.pdf", width = 7.5, height = 4.2, units='in', dpi=600)

############################################################### Figure A5b

p2_1 <- marginaleffect_logit2(ModelResults =  model_onset_relative_TEK, n.sim = 1000, 
                              varname = c("best_low_tlineq", "worst_low_tlineq",
                                          "median_low_tlineq","nearst_low_tlineq"),
                              data = data_tek,
                              dv = "relative_poverty",
                              val1 = 0, val2 =1, clusterid = "gwgroupid") 


p2_2 <- marginaleffect_logit2(ModelResults =  model_onset_relative_TEK, n.sim = 1000, 
                              varname = c("best_high_tlineq", "worst_high_tlineq",
                                          "median_high_tlineq","nearst_high_tlineq"),
                              data = data_tek,
                              dv = "relative_wealth",
                              val1 = 0, val2 =1, clusterid = "gwgroupid") 


df_relative_onset <- bind_rows(p2_1, p2_2)



p = ggplot(df_relative_onset , aes(colour = dv)) + 
       geom_hline(yintercept = 0, lty = 2) +
       geom_linerange(aes(x = id, ymin = low,
                          ymax = high),
                      lwd = 1.5, position = position_dodge(width = 1/2)) +
       geom_pointrange(aes(x = id, y =median, ymin =low,
                           ymax =high, shape = dv),
                       lwd = 1.5, position = position_dodge(width = 1/2),
                       fill = "WHITE")+
       geom_text(aes(x = id +0.1 ,
                     y = median,
                     label = format(round(median, 2))),position = position_dodge(width = 1/2)) +
       theme_bw() + xlab("") + ylab("First Difference in Predicted Prob.") + 
       theme(legend.position="bottom",
             legend.title=element_blank(),
             legend.text = element_text(colour="blue", size=12),
             axis.text= element_text(size=12),
             text = element_text(size=12),
             plot.title = element_text(hjust = .5, size = 12))+
       ggtitle("Conflict Onset(More than One TEK)")+
       scale_colour_discrete(labels = rev(c("Relative Wealth","Relative poverty")))+
       scale_shape_discrete(labels = rev(c("Relative Wealth","Relative poverty")) )+
       guides(color = guide_legend(nrow = 1))+
       scale_x_continuous(breaks = 1:4,labels =parse(text =rev(c( 
              "Nearest~TEK",
              "Median~TEK",
              "Worst~TEK",
              "Best~TEK")))) 

ggsave("Final_JPR/Figures/appendix_figures/fig_A5-b.pdf", width = 7.5, height = 4.2, units='in', dpi=600)



## Figure A6

###################################################################################
## drop Ira
data_iraq <- data %>% dplyr::filter(countries_gwid != 645)
iraq <- data %>% dplyr::filter(countries_gwid == 645)


dv <- 'onset_do_flag' 
## controls
ctrs <- c('lineq_total', 'status_excl','lsize','family_warhistdummy','family_downgraded2', 
          'ctygdppc_log', 'ctypop_log', 'peaceyears','peaceyears_sq','peaceyears_cub')

ivs_1 <- c('best_tlineq_total', 'tekbest_status_excl')
ivs_2 <- c('worst_tlineq_total', 'tekworst_status_excl')
ivs_3 <- c("median_tlineq_total",'tekmedian_status_excl')
ivs_4 <- c('nearst_tlineq_total','teknearst_status_excl')

f_do_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))

#
ivs <- c('best_tlineq_total','worst_tlineq_total','median_tlineq_total','nearst_tlineq_total')


## using log of transnational inequality
model_onset_iraq <- list()
for (i in 1:length(ivs)){
       model_onset_iraq [[i]] <-glm(eval(parse(text = paste0('f_do_', i))), 
                                    data = data_iraq, family = binomial(link = "logit"))
}


## dependent variable
dv <- 'incidence_flag' 
## controls
ctrs <- c('lineq_total', 'status_excl','lsize','family_warhistdummy','family_downgraded2', 
          'ctygdppc_log', 'ctypop_log', 'peaceyears','peaceyears_sq','peaceyears_cub')

ivs_1 <- c('best_tlineq_total', 'tekbest_status_excl')
ivs_2 <- c('worst_tlineq_total', 'tekworst_status_excl')
ivs_3 <- c("median_tlineq_total",'tekmedian_status_excl')
ivs_4 <- c('nearst_tlineq_total','teknearst_status_excl')

f_do_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))

#
ivs <- c('best_tlineq_total','worst_tlineq_total','median_tlineq_total','nearst_tlineq_total')


## using log of transnational inequality
model_incidence_iraq <- list()
for (i in 1:length(ivs)){
       model_incidence_iraq[[i]] <-glm(eval(parse(text = paste0('f_do_', i))), 
                                       data = data_iraq, family = binomial(link = "logit"))
}



p1_1 <- marginaleffect_logit2(ModelResults =   model_onset_iraq, n.sim = 1000, 
                              varname = c("best_tlineq_total", "worst_tlineq_total",
                                          "median_tlineq_total","nearst_tlineq_total"),
                              data = data_iraq,
                              dv = "onset",
                              val1 = 0, val2 =1, clusterid = "gwgroupid") 


p1_2 <- marginaleffect_logit2(ModelResults =  model_incidence_iraq , n.sim = 1000, 
                              varname = c("best_tlineq_total", "worst_tlineq_total",
                                          "median_tlineq_total","nearst_tlineq_total"),
                              data = data_iraq,
                              dv = "incidence",
                              val1 = 0, val2 =1, clusterid = "gwgroupid") 


df_iraq <- bind_rows(p1_1, p1_2)


## Figure A6a

p = ggplot(df_iraq , aes(colour = dv)) + 
       geom_hline(yintercept = 0, lty = 2) +
       geom_linerange(aes(x = id, ymin = low,
                          ymax = high),
                      lwd = 1.5, position = position_dodge(width = 1/2)) +
       geom_pointrange(aes(x = id, y =median, ymin =low,
                           ymax =high, shape = dv),
                       lwd = 1.5, position = position_dodge(width = 1/2),
                       fill = "WHITE")+
       geom_text(aes(x = id +0.1 ,
                     y = median,
                     label = format(round(median, 2))),position = position_dodge(width = 1/2)) +
       theme_bw() + xlab("") + ylab("First Difference in Predicted Prob.") + 
       theme(legend.position="bottom",
             legend.title=element_blank(),
             legend.text = element_text(colour="blue", size=12),
             axis.text= element_text(size=12),
             text = element_text(size=12),
             plot.title = element_text(hjust = .5, size = 12))+
       ggtitle("Logistic Models (Excluding Iraq)")+
       scale_colour_discrete(labels = rev(c("Conflict Onset","Conflict Incidence")))+
       scale_shape_discrete(labels = rev(c("Conflict Onset","Conflict Incidence")) )+
       guides(color = guide_legend(nrow = 1))+
       scale_x_continuous(breaks = 1:4,labels =parse(text =rev(c( 
              "Nearest~TEK",
              "Median~TEK",
              "Worst~TEK",
              "Best~TEK")))) 

ggsave("Final_JPR/Figures/appendix_figures/fig_A6-a.pdf", width = 7.5, height = 4.2, units='in', dpi=600)


## Figure A6-c
# relative wealth exlcuding Iraq
dv <- 'onset_do_flag' 
## controls
ctrs <- c('lineq_total', 'status_excl','lsize','family_warhistdummy','family_downgraded2', 
          'ctygdppc_log', 'ctypop_log','peaceyears','peaceyears_sq','peaceyears_cub')

ivs_1 <- c('best_low_tlineq', 'best_high_tlineq', 'tekbest_status_excl')
ivs_2 <- c('worst_low_tlineq', 'worst_high_tlineq', 'tekworst_status_excl')
ivs_3 <- c('median_low_tlineq', 'median_high_tlineq','tekmedian_status_excl')
ivs_4 <- c('nearst_low_tlineq', 'nearst_high_tlineq','teknearst_status_excl')


f_do_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))

#
ivs <- c('best','worst','median','nearst')


## using log of transnational inequality
model_onset_relative_iraq <- list()
for (i in 1:length(ivs)){
       model_onset_relative_iraq[[i]] <- glm(eval(parse(text = paste0('f_do_', i))), 
                                             data = data_iraq, family = binomial(link = "logit"))
}


## dependent variable
dv <- 'incidence_flag' 
## controls
ctrs <- c('lineq_total', 'status_excl','lsize','family_warhistdummy','family_downgraded2', 
          'ctygdppc_log', 'ctypop_log','peaceyears','peaceyears_sq','peaceyears_cub')

ivs_1 <- c('best_low_tlineq', 'best_high_tlineq', 'tekbest_status_excl')
ivs_2 <- c('worst_low_tlineq', 'worst_high_tlineq', 'tekworst_status_excl')
ivs_3 <- c('median_low_tlineq', 'median_high_tlineq','tekmedian_status_excl')
ivs_4 <- c('nearst_low_tlineq', 'nearst_high_tlineq','teknearst_status_excl')


f_do_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_do_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))

#
ivs <- c('best','worst','median','nearst')


## using log of transnational inequality
model_incidence_relative_iraq <- list()
for (i in 1:length(ivs)){
       model_incidence_relative_iraq[[i]] <- glm(eval(parse(text = paste0('f_do_', i))), 
                                                 data = data_iraq, family = binomial(link = "logit"))
}


p1_1 <- marginaleffect_logit2(ModelResults =  model_incidence_relative_iraq, n.sim = 1000, 
                              varname = c("best_low_tlineq", "worst_low_tlineq",
                                          "median_low_tlineq","nearst_low_tlineq"),
                              data = data_iraq,
                              dv = "relative_poverty",
                              val1 = 0, val2 =1, clusterid = "gwgroupid") 


p1_2 <- marginaleffect_logit2(ModelResults =  model_incidence_relative_iraq, n.sim = 1000, 
                              varname = c("best_high_tlineq", "worst_high_tlineq",
                                          "median_high_tlineq","nearst_high_tlineq"),
                              data = data_iraq,
                              dv = "relative_wealth",
                              val1 = 0, val2 =1, clusterid = "gwgroupid") 


df_relative_incidence_iraq <- bind_rows(p1_1, p1_2)



p = ggplot(df_relative_incidence_iraq , aes(colour = dv)) + 
       geom_hline(yintercept = 0, lty = 2) +
       geom_linerange(aes(x = id, ymin = low,
                          ymax = high),
                      lwd = 1.5, position = position_dodge(width = 1/2)) +
       geom_pointrange(aes(x = id, y =median, ymin =low,
                           ymax =high, shape = dv),
                       lwd = 1.5, position = position_dodge(width = 1/2),
                       fill = "WHITE")+
       geom_text(aes(x = id +0.1 ,
                     y = median,
                     label = format(round(median, 2))),position = position_dodge(width = 1/2)) +
       theme_bw() + xlab("") + ylab("First Difference in Predicted Prob.") + 
       theme(legend.position="bottom",
             legend.title=element_blank(),
             legend.text = element_text(colour="blue", size=12),
             axis.text= element_text(size=12),
             text = element_text(size=12),
             plot.title = element_text(hjust = .5, size = 12))+
       ggtitle("Conflict Incidence(Excluding Iraq)")+
       scale_colour_discrete(labels = rev(c("Relative Wealth","Relative poverty")))+
       scale_shape_discrete(labels = rev(c("Relative Wealth","Relative poverty")) )+
       guides(color = guide_legend(nrow = 1))+
       scale_x_continuous(breaks = 1:4,labels =parse(text =rev(c( 
              "Nearest~TEK",
              "Median~TEK",
              "Worst~TEK",
              "Best~TEK")))) 

ggsave("Final_JPR/Figures/appendix_figures/fig_A6-c.pdf", width = 7.5, height = 4.2, units='in', dpi=600)

### Figure A6-b

p2_1 <- marginaleffect_logit2(ModelResults =  model_onset_relative_iraq, n.sim = 1000, 
                              varname = c("best_low_tlineq", "worst_low_tlineq",
                                          "median_low_tlineq","nearst_low_tlineq"),
                              data = data_iraq,
                              dv = "relative_poverty",
                              val1 = 0, val2 =1, clusterid = "gwgroupid") 


p2_2 <- marginaleffect_logit2(ModelResults =  model_onset_relative_iraq, n.sim = 1000, 
                              varname = c("best_high_tlineq", "worst_high_tlineq",
                                          "median_high_tlineq","nearst_high_tlineq"),
                              data = data_iraq,
                              dv = "relative_wealth",
                              val1 = 0, val2 =1, clusterid = "gwgroupid") 


df_relative_onset <- bind_rows(p2_1, p2_2)



p = ggplot(df_relative_onset , aes(colour = dv)) + 
       geom_hline(yintercept = 0, lty = 2) +
       geom_linerange(aes(x = id, ymin = low,
                          ymax = high),
                      lwd = 1.5, position = position_dodge(width = 1/2)) +
       geom_pointrange(aes(x = id, y =median, ymin =low,
                           ymax =high, shape = dv),
                       lwd = 1.5, position = position_dodge(width = 1/2),
                       fill = "WHITE")+
       geom_text(aes(x = id +0.1 ,
                     y = median,
                     label = format(round(median, 2))),position = position_dodge(width = 1/2)) +
       theme_bw() + xlab("") + ylab("First Difference in Predicted Prob.") + 
       theme(legend.position="bottom",
             legend.title=element_blank(),
             legend.text = element_text(colour="blue", size=12),
             axis.text= element_text(size=12),
             text = element_text(size=12),
             plot.title = element_text(hjust = .5, size = 12))+
       ggtitle("Conflict Onset(Excluding Iraq)")+
       scale_colour_discrete(labels = rev(c("Relative Wealth","Relative poverty")))+
       scale_shape_discrete(labels = rev(c("Relative Wealth","Relative poverty")) )+
       guides(color = guide_legend(nrow = 1))+
       scale_x_continuous(breaks = 1:4,labels =parse(text =rev(c( 
              "Nearest~TEK",
              "Median~TEK",
              "Worst~TEK",
              "Best~TEK")))) 
ggsave("Final_JPR/Figures/appendix_figures/fig_A6-b.pdf", width = 7.5, height = 4.2, units='in', dpi=600)


